library(readxl)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
feature_id <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "feature_id")
features <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "features")
metadata <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "metadata")
library(dplyr)
library(tibble)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
diff_all <- features %>% column_to_rownames("features") %>% t() %>% as.data.frame() %>% rownames_to_column(var="sample_ID") %>% merge(metadata, by="sample_ID")
matrix_diff <- diff_all %>% remove_rownames() %>% column_to_rownames(var="sample_ID")
matrix_diff=as.matrix(matrix_diff)
class(matrix_diff) <- "numeric"
Warning: NAs introduced by coercion
matrix_diff[is.na(matrix_diff)] = 0
#ANNOTATION COLUMNS
anno_row <- diff_all %>%
dplyr::select(Treatment,Time_tx)
rownames(anno_row)= rownames(as.data.frame(matrix_diff))
anno_row <- as.data.frame(anno_row)
anno_color <- list(Treatment = c(Antibiotic = "#DF8F44FF", Control = "#374E55FF"),
Time_tx = c(`Week 1` = "#EE0000FF", `Day -1` = "#3B4992FF", `Week 9` = "#008B45FF"))
library(pheatmap)
library(viridis) #color pallet, it's optional
library(RColorBrewer)
#HEATMAP RA USING ln color pallet RBrewer
heatmap_metIDs <- pheatmap(
#mat = t(matrix_diff),
mat = log10(t(matrix_diff+0.0000001)),
border_color = NA,
show_colnames = T,
show_rownames = T,
angle_col = 90,
drop_levels = TRUE,
fontsize_col = 4,
fontsize_row = 5,
fontsize = 14,
color = brewer.pal(9,"Reds"),
number_color = NA,
annotation_col = anno_row,
annotation_colors = anno_color,
annotation_names_col = T,
annotation_names_row = T,
cluster_cols = T,
cluster_rows = T,
clustering_method = "ward.D",
gaps_row = FALSE
)
Warning: NaNs produced
heatmap_metIDs
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
matrix_features <- read_excel("metabolomics_microbiome_fc.xlsx", sheet = "features")%>% remove_rownames %>% column_to_rownames(var="features")
matrix_features=as.matrix(matrix_features)
class(matrix_diff) <- "numeric"
matrix_diff[is.na(matrix_diff)] = 0
hc <- hclust(dist(matrix_features), "ward.D")
plot(hc, cex =.5)
hc_groups <- as.data.frame(cutree(hc, k = 30)) %>% tibble::rownames_to_column("features") %>% rename(hc_group = 2)
hc_groups %>% count(hc_group)
#Number of clusters
library("NbClust")
nb <- NbClust(diss = dist(matrix_features), distance = NULL, min.nc = 2,
max.nc = 50, method = "ward.D",index = "silhouette")
Only frey, mcclain, cindex, sihouette and dunn can be computed. To compute the other indices, data matrix is needed
library(janitor)
library(tidyr)
diff <- as.data.frame(matrix_features) %>% rownames_to_column("features")
df_hc_merge <- merge(diff, hc_groups, by = "features")
df_hc <- df_hc_merge %>% gather(sample_ID, abundance, -c(features,hc_group))
#metadata <- metadata %>% rename(sample_ID = ID)
df_hc_meta <- merge(df_hc, metadata, by="sample_ID")
df_hc_meta$abundance <- as.numeric(as.character(df_hc_meta$abundance))
df_hc_meta_type <- merge(df_hc_meta, feature_id, by="features")
hc2 <- df_hc_meta_type %>%
filter(hc_group==2, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,15)
hc2
hc4 <- df_hc_meta_type %>%
filter(hc_group==4, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,5)
hc4
hc5 <- df_hc_meta_type %>%
filter(hc_group==5, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc5
hc6 <- df_hc_meta_type %>%
filter(hc_group==6, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,5)
hc6
hc7 <- df_hc_meta_type %>%
filter(hc_group==7, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc7
hc8 <- df_hc_meta_type %>%
filter(hc_group==8, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc8
hc9 <- df_hc_meta_type %>%
filter(hc_group==9, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,10)
hc9
hc10 <- df_hc_meta_type %>%
filter(hc_group==10, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc10
hc11 <- df_hc_meta_type %>%
filter(hc_group==11, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+
theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc11
hc12 <- df_hc_meta_type %>%
filter(hc_group==12, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+
theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc12
hc14 <- df_hc_meta_type %>%
filter(hc_group==14, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx",panel.labs = list(Time_tx=c("Lactation","Dry-off","Fresh cows"))) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#DF8F44FF","#B24745FF")) + scale_color_manual(values=c("#DF8F44FF","#B24745FF"))
hc14
hc15 <- df_hc_meta_type %>%
filter(hc_group==15, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+
theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,10)
hc15
hc28 <- df_hc_meta_type %>%
filter(hc_group==14, Time_tx%in%c("Day -1","Week 1","Week 9")) %>%
ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+
theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc28
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")
#ggsave(plot=hc1,"hc1_multiomics.png",width = 16, height = 8)
ggsave(plot=hc2,"hc2_multiomics.png",width = 30, height = 20)
#ggsave(plot=hc3,"hc3_multiomics.png",width = 20, height = 20)
#ggsave(plot=hc4,"hc4_multiomics.png",width = 16, height = 8)
ggsave(plot=hc5,"hc5_multiomics.png",width = 30, height = 20)
#ggsave(plot=hc14,"hc14_multiomics.png",width = 16, height = 20)
#ggsave(plot=hc15,"hc15_multiomics.png",width = 16, height = 20)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")
hc_bars <- ggarrange(hc1,hc14,hc3,labels = c("A","B","C"),nrow = 3, ncol = 1)
ggsave(plot=hc_bars,"hc_bars_fc_multiomics.png",width = 15, height = 15)